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ABSTRACT 

This  report  presents  the  development  and  automated  numerical  implementation  of  an 
iterative  gradientless  optimisation  method  for  the  analysis  of  problems  relating  to  life 
extension  of  aircraft  components.  The  method  has  been  implemented  to  interface  with 
the  finite  element  code  PAFEC,  which  does  not  normally  have  an  optimisation 
capability.  The  key  feature  of  the  approach  is  to  achieve  constant  boxmdary  stresses,  in 
regions  of  interest,  by  moving  nodes  on  the  stress  concentrator  boxmdary  by  an 
amotmt  dependent  on  the  sign  and  magnitude  of  the  local  hoop  stress  obtained  from  a 
previous  iteration  of  a  standard  finite  element  analysis.  The  results  of  example 
problems  are  presented  which  include  the  optimisation  of  hole  shapes  in  flat  plates 
and  the  optimisation  of  the  design  of  bonded  reinforcements  with  a  focus  on 
minimising  adhesive  stress  while  maintaining  the  effectiveness  of  the  reinforcement. 
In  aU  cases  significant  stress  reductions  were  achieved  by  way  of  the  local  shape 
changes.  The  method  presented  is  considered  a  simple  robust  complementary  method 
to  the  use  of  commercially  available  gradient  based  finite  element  optimisation 
software.  It  is  also  considered  suitable  for  use  with  typical  standard  commercial  finite 
element  packages  other  than  PAFEC. 
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Executive  Summary 

Two  common  fatigue  life  extension  practices  which  seek  to  favourably  alter  the  local 
stress  distribution  near  a  stress  concentrator  in  a  RAAF  aircraft  component  are:  (i) 
modification  of  the  existing  local  geometry  by  removing  material,  and  (ii)  addition  of 
new  material  to  the  local  region;  such  as  the  adhesive  bonding  of  a  composite  patch.  In 
both  these  cases  an  improved  geometry  (i.e.  shape)  is  sought  that  minimises  the 
magnitude  of  peak  terrsile  stresses.  However,  the  design  of  re-work  shapes  or  bonded 
patches  for  RAAF  aircraft  components  has  typically  been  based  on  intuitive  trial  and 
error  (or  simplified  analytical  formulations).  Subsequently  the  effectiveness  of  such 
design  has  usually  been  quantified  using  standard  finite  element  analyses  methods.  It 
is  believed  that  the  use  of  optimisation  procedures  would  significantly  improve  the 
effectiveness  of  such  shape  modification  designs.  However,  finite  element  based 
optimisation  procedures  have  not  previously  been  used  at  AMRL  in  this  context. 

This  report  presents  the  development  and  automated  numerical  implementation  of  an 
iterative  gradientless  optimisation  method  for  the  analysis  of  problems  relating  to  life 
extension  of  aircraft  components.  The  method  has  been  implemented  to  interface  with 
the  finite  element  code  PAFEC,  which  does  not  normally  have  an  optimisation 
capability.  The  key  featiue  of  the  approach  is  to  achieve  constant  boundary  stresses,  in 
regions  of  interest,  by  moving  nodes  on  the  stress  concentrator  boundary  by  an 
amount  dependent  on  the  sign  and  magnitude  of  the  local  hoop  stress  obtained  from  a 
previous  iteration  of  a  standard  finite  element  analysis.  The  results  of  example 
problems  are  presented  which  include  the  optimisation  of  hole  shapes  in  flat  plates 
and  the  optimisation  of  the  design  of  bonded  reinforcements  with  a  focus  on 
minimising  adhesive  stress  while  maintaining  the  effectiveness  of  the  reinforcement. 
In  all  cases  significant  stress  reductions  were  achieved  by  way  of  the  local  shape 
changes.  The  method  presented  is  considered  a  simple  robust  complementary  method 
to  the  use  of  commercially  available  gradient  based  finite  element  optimisation 
software. 

The  procedure  developed  and  implemented  in  this  report  has  given  AED  the 
capability  of  determining  optimal  shapes  for  extending  the  fatigue  hfe  of  critical 
regions  in  RAAF  aircraft  structural  components. 
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Notation 


stress  in  tangential  (hoop)  direction 
contour  of  nodes  on  bovmdary 
shear  stress  in  x,  y  co-ordinate  system 
Poisson's  Ratio 

movement  of  node  in  normal  direction 

elastic  modulus 

gigapascal 

lengA  of  tapered  region  of  patch 

millimetre 

megapascal 

node  number 

total  number  of  nodes  on  boundary 

local  cartesion  co-ordinates  where  n  is  normal  to  contour 

characteristic  length 

step  size  scahng  factor 

patch  thickness  in  tapered  region  as  a  fxmction  of  x 

global  cartesian  co-ordinates 

objective  fimction  to  monitor  convergence 


Subscripts 


node  number 

total  number  of  nodes  on  boimdary 
refers  to  initial  patch  thickness 
threshold  value 
average  value 


Superscripts 

value  after  smoothing 
maximum  value 
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AMRL 

FORTRAN 

FFVH 

FFVH13 

MSC/NASTRAN 

PAFEC 

PAFBLOCK 

RAAF 


Aeronautical  and  Maritime  Research  Laboratory 
computer  language  for  FORmula  TRANslation 
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1.  Introduction 


AMRL  has  a  continuing  focus  on  implementing  and  developing  methods  suitable  for 
the  fatigue  life  extension  of  RAAF  aircraft  structural  components.  These  methods  are 
either  directed  at  delaying  the  onset  of  crack  initiation,  and/  or  reducing  the  rate  of 
crack  growth  at  critical  stress  concentrating  locations.  These  critical  locations  are  often 
on  free  boundaries  of  structural  components,  and  have  a  local  geometry  such  that  they 
are  subjected  to  elevated  cyclic  tensile  stresses  when  the  airframe  undergoes  service 
loading.  Two  common  life  extension  practices  which  seek  to  alter,  in  a  favourable 
manner,  the  local  stress  distribution  near  a  stress  concentrator  are:  (i)  modification  of 
the  existing  local  geometry  by  removing  material,  and  (ii)  addition  of  new  material  to 
the  local  region;  for  example  by  using  an  adhesively  bonded  reinforcement  or  an 
interference  fitted  component.  In  both  these  cases  an  improved  geometry  (i.e.  shape) 
that  minimises  the  peak  magnitude  of  tensile  stresses  at  a  stress  concentrating  local 
region  is  desirable. 

Finite  element  analysis  provides  a  powerful  tool  for  determining  the  stresses  for  a 
given  geometry,  and  potentially  for  determining  an  optimised  geometry  to  minimise 
the  stress  concentration.  However,  while  finite  element  analyses  have  been  extensively 
used  at  AMRL,  it  is  important  to  note  that  finite  element  based  optimisation 
procedures  have  not  previously  been  applied  at  AMRL.  In  the  past,  rework  shapes 
have  typically  been  developed  by  using  intuitive  trial  and  error,  and  then  their 
effectiveness  quantified  using  standard  finite  element  analyses.  Sunilarly,  in  the  case  of 
bonded  repairs,  the  initial  design  has  typically  been  based  on  simplified  analytical 
formulations,  and  then  finalised  by  imdertaking  standard  finite  element  analyses. 
Often  the  design  of  the  bonded  repairs  has  been  confined  to  rectangular  patches  of 
constant  thickness  with  linear  tapering  at  the  ends. 

A  number  of  finite  element  based  approaches  have  been  given  in  the  literature  for 
shape  optimisation  to  minimise  stress  concentrations.  In  aU  these  approaches  the  initial 
geometry  is  modified  between  iterations  of  a  number  of  finite  element  analyses,  until 
convergence  of  the  peak  stress  to  a  suitable  value  has  been  achieved.  A  number  of 
procedures  can  be  employed  for  iterating  the  shape,  and  these  include  sensitivity 
analysis  methods  [1-5],  gradientless  methods  [6-10]  and  layout  methods  [11-13]. 

In  this  paper  we  develop  a  simple,  efficient  gradientless  computational  method  tiiat 
can  be  used  with  any  standard  commercial  finite  element  package  for  minimising 
stress  concentrations.  The  driving  force  for  the  method  is  to  achieve  constant  boimdary 
stresses  by  correctly  moving  boimdary  nodal  positions,  in  this  sense  it  has  most  in 
common  with  the  work  of  Schnack  [6,7\  and  Mattheck  et  al  [8-10].  However,  there  are 
significant  differences,  primarily  in  the  algorithm  for  determining  new  nodal  position 
between  iterations.  In  Section  2,  a  brief  discussion  of  various  finite  element  based 
shape  optimisation  approaches,  including  sensitivity,  gradientiess  and  layout 
me^ods,  are  given.  In  Section  3  a  moving  boundary  method  is  postulated  and 
presented  for  gradientless  shape  optimisation,  using  only  the  results  of  standard  finite 
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element  methods,  in  an  iterative  manner.  In  this  approach  nodes  on  the  stress 
concentrator  boundary  are  moved  by  an  amount  dependent  on  the  sign  and 
magnitude  of  the  local  hoop  stress,  obtained  from  a  previous  iteration  of  a  standard 
finite  element  analysis.  The  results  obtained  from  application  of  the  method  to  sample 
problems  involving  stress  concentrators  in  two  dimensional  plates  are  presented  in 
Section  4.  Subsequently  in  Section  5,  numerical  examples  relating  to  bonded  repairs  in 
two  dimensions  are  given.  In  the  appendix  software  for  the  automated  implementation 
of  the  procedure  is  presented  for  use  with  the  PAFEC  finite  element  package. 


2.  Comments  on  Finite  Element  Based  Shape 
Optimisation  Methods 


2.1  Sensitivity  analysis  methods 

Approaches  based  on  sensitivity  analysis  require  stress  gradient  information  to  be 
determined,  so  that  the  sensitivity  of  ^e  objective  fxmction  to  boimdary  movements 
(i.e.  shape  change)  at  the  stress  concentrator  can  be  quantified  [1-5].  Most  industry- 
standard  finite  element  analysis  software  codes  which  have  an  optimisation  capability 
use  this  approach.  However,  many  commercial  finite  element  codes  have  no 
optimisation  capability.  While  sensitivity  analysis  is  a  powerful  general  method  for 
optimisation,  it  can  also  be  relatively  demanding  on  computer  time.  Also  sigrtificant  is 
the  fact  that  typically  a  great  deal  of  effort  is  required  (compared  to  standard  finite 
element  analyses)  to  generate  a  suitable  mitial  model  specification  to  obtain  a 
successful  result,  with  mesh  distortion  and  appropriate  specification  of  the  objective 
function  being  particular  problems.  In  the  open  literature  a  number  of  papers  deal 
with  the  development  of  this  method  and  apply  it  to  relatively  simple  geometries.  A 
particular  focus  has  been  on  developing  methods  to  reduce  computational  time. 
Surprisingly,  to  date  very  few  examples  are  available  in  the  literature  discussing 
applications  to  more  complex  geometries. 


2.2  Gradientless  analysis  methods 

In  gradientless  analysis  approaches  stress  derivatives  are  not  used  to  determine 
optimal  geometries,  hence  they  are  typically  much  simpler  to  implement,  or  use  with 
standard  (i.e.  non-optimisation)  finite  element  codes.  The  driving  force  for  these 
methods  is  to  determine  a  bormdary  shape,  such  that  there  is  a  constant  or  near 
constant  stress  distribution  around  the  boimdary.  Previous  work  [6-10,14]  has  shown 
that  such  a  situation  results  in  a  minimised  peak  stress  aroimd  the  boimdary,  except 
for  the  presence  of  structural  details  that  restrict  shape  changes.  Schnack  et  al  [6,7]  and 
Mattheck  et  al.  [8-10]  in  particular  have  used  finite  element  based  methods  to 
determinal  optimal  shape  changes.  However,  the  logic  that  these  two  authors 
respectively  use  to  determine  the  new  shape  is  completely  different. 
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The  work  of  Mattheck  et  al.  is  based  on  the  interesting  observations  of  shape 
optimisation  in  living  organisms  such  as  in  trees  and  animal  bones.  As  part  of  the 
evolutionary  process  many  such  'living  structures'  appear  to  have  been  able  to  add 
material  in  regions  of  high  stress  and  reduce  material  in  regions  of  low  stress  to  bring 
about  an  optimal  shape  that  produces  a  constant  Von  hlises  stress  distribution  on  the 
free  surface.  To  mimic  this  process  for  engineering  structures,  he  uses  an  iterative 
approach  as  follows.  First,  undertake  a  standard  finite  element  nm  to  determine 
stresses  on  the  boundary.  Second,  do  a  subsequent  analysis  in  which  an  adhesive  layer 
is  added  to  the  boundary,  and  is  allowed  to  swell  or  contract  in  proportion  due  to  the 
previously  determined  boundary  stresses  (which  have  been  converted  to  an  equivalent 
temperature  distribution).  The  amoimt  of  swelling  or  contraction  of  the  adhesive  layer 
determines  where  to  add  material  or  remove  material  in  the  next  fmite  element 
analysis,  which  is  where  the  next  iteration  occurs. 

In  the  work  of  Schnack  et  al,  an  extension  to  the  fade  away  hypothesis  of  Neuber,  for 
minimising  peak  notch  stresses  is  used.  There  are  three  key  ideas  in  Schnack's  work. 
Firstly,  that  a  change  in  the  curvature  at  a  point  on  a  notch  surface  has  a  significant 
effect  on  the  magnitude  of  the  tangential  stress  at  that  point.  Secondly,  the  higher  the 
peak  stress  at  a  point  on  the  notch,  the  more  the  stresses  will  fade  away  as  the  distance 
from  the  high  stress  zone  increases.  Thirdly,  the  magnitude  of  the  maximum  stress  can 
be  reduced  by  increasing  the  minimum  stress  on  the  notch  surface.  Schnack's  finite 
element  implementation  starts  by  completing  a  standard  analysis  and  determining 
tangential  stresses  at  location  around  the  hole  boimdary.  Then  the  bormdary  nodes 
near  the  location  of  maximum  tangential  stress  are  moved  to  increase  the  local  radius 
of  curvature  at  that  point.  Similarly  the  nodes  near  the  location  of  minimum  tangential 
stress  are  also  moved  to  decrease  the  local  curvature.  At  the  same  time  aU  other  nodes 
around  the  boundary  are  moved  to  achieve  a  smooth  transition  of  curvature  between 
these  two  extreme  points.  With  these  new  nodal  positions,  another  finite  element 
analysis  is  then  rmdertaken  and  the  process  repeated  until  convergence  is  achieved. 


2.3  Layout  analysis  methods 

In  typical  layout  optimisation  [11-13],  a  highly  discretised  finite  element  mesh  is  used, 
of  uniform  element  size,  for  a  given  design  space.  Then,  subject  to  fixed  geometry 
constraints,  hghtly  stressed  elements  from  a  finite  element  model  are  removed  or  as  a 
minor  variation,  the  stiffness  of  lightly  stressed  elements  is  successively  reduced.  This 
approach  then  leads  to  the  generation  of  many  voids,  so  that  in  the  limit  open  frame, 
truss-like  structures  can  be  created  which  are  approximately  uniformly  stressed. 
However,  if  desired,  the  generation  of  multiple  holes  can  be  suppressed  [13],  so  that 
only  the  geometry  of  one  hole  is  modified  (i.e.  by  removing  only  elements  only  that  are 
at  the  hole  boxmdary).  Typically,  although  the  general  hole  shape  will  be  correct,  the 
boundary  will  'jagged'.  Therefore  another  finite  element  analysis  is  then  needed  (using 
a  fitted  smooth  boimdary)  to  determine  the  appropriate  stress  concentrations.  Hence 
these  methods,  although  generally  not  able  to  generate  smooth  boundaries,  appear 
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effective  in  determining  an  approximate  shape,  which  can  be  used  if  desired  as  a 
starting  geometry  for  analysis  by  other  methods. 


3.  Moving  Boundary  Method 

3.1  General  description 

Here  we  aim  to  develop  a  simple,  efficient  gradientless  computational  method  that  can 
be  used  with  any  standard  commercial  finite  element  package  for  minimising  stress 
concentrations.  The  driving  force  for  the  metiiod  is  to  achieve  constant  boundary  stress 
by  correctly  moving  boundary  nodal  positions,  in  this  sense  the  present  method  has 
most  in  common  with  the  works  of  Schnack  [6,7]  and  Mattheck  et  al  [8-10]  and  their  co¬ 
workers.  However,  there  are  significant  differences  in  the  present  approach,  primarily 
in  the  algorithm  for  determining  new  nodal  positions  between  iterations,  as  detailed  in 
the  following  sections. 

Consider  Figure  1  which  shows  a  hole  (stress  concentrator)  with  boundary,  T  on 
which  there  are  a  number  of  nodes  i  =  l,k,  under  an  arbitrary  remote  loading.  The 
tangential  stress  will  vary  around  the  boimdary  due  to  remote  loading  as  shown 
schematically  in  Figure  2.  Let  us  use  the  logic  of  Mattheck  (and  others)  and  add 
material  where  stresses  are  high  and  remove  it  where  stresses  are  low.  The  question  is 
then  what  should  the  relative  amounts  of  material  to  be  added  or  removed  be.  A 
simple  postulate  would  be  to  require  the  amount  of  material  added  or  removed  at  any 
point  on  the  boimdary  to  be  directly  proportional  to  the  difference  between  the  local 
tangential  stress  and  a  suitable  reference  value,  and  this  process  can  then  be  repeated 
iteratively  until  a  reduced  stress  concentration  has  been  achieved.  Hence  the  amoimt 
to  move  a  given  node  i  on  the  boundary  (in  the  normal  direction)  can  be  written  as, 

d.  =  (1) 

^th 

where  positive  d,-  indicates  material  addition,  cr,-  is  the  tangential  stress  at  node  i  on 
the  boundary,  cr,^  is  a  threshold  boundary  hoop  stress,  r  is  an  arbitrary  characteristic 
length  such  as  initial  hole  radius,  and  s  is  an  arbitrary  step  size  scaling  factor  typically 
in  the  range  of  0.1-0.4. 

A  number  of  interesting  features  immediately  become  apparent,  particularly  in 
relation  to  the  selection  of  the  parameter  cr,^,  as  follows;  (i)  if  we  select  cr less  than 
the  maximum,  but  more  than  the  minimum  stress  determined  from  the  previous 
iteration,  material  will  be  both  added  and  removed  at  various  locations  around  the 
boundary,  (ii)  if  we  select  a  the  same  or  greater  than  the  maximum  stress 
determined  from  the  previous  iteration,  the  boundary  will  be  moved  to  remove 
material  only,  and  (iii)  if  we  select  to  be  the  same  or  less  than  the  minimum  stress 
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determined  from  the  previous  iteration,  the  botmdary  will  be  moved  to  add  material 
only.  If  we  are  reworking  a  stress  concentrator  in  an  existing  component,  we  must  use 
option  (ii)  or  option  (i)  with  a  further  constraint  on  botmdary  movements  so  that 
material  is  not  added  (i.e.  d,-  <  0 ). 

This  method  is  suitable  for  implementation  with  a  standard  finite  element  code  that 
allows  for  movement  of  nodes  after  the  completion  of  an  analysis.  The  procedure  is 
shown  in  the  simple  flow  chart  given  in  Figure  3.  For  such  an  iterative  procedure  an 
appropriate  objective  function,  Y,  to  monitor  solution  convergence  can  be  written  as 


Y  = 


k-i 


(2) 


Clearly,  as  7  approaches  zero,  the  stresses  become  more  uniform  aroimd  the 
boxmdary.  It  should  be  noted  that  a  completely  uniform  distribution  of  stresses  around 
the  boundary  can  only  be  achieved  for  a  restricted  set  of  problems.  The  approach  was 
implemented  with  the  finite  element  code  PAFEC,  and  this  is  described  in  the 
Appendix  A.  It  should  be  noted  that  in  the  implementation  with  PAFEC,  only  the 
nodal  movements  at  the  stress  concentrator  botmdary  need  to  be  specified.  Internal 
nodes  are  automatically  moved  in  a  proportionate  manner  by  making  use  of  an 
element  generation  feature  named  PAFBLOCKs.  PAFBLOCKs  are  typically  four  sided 
subregioits  in  the  finite  element  model  consisting  of  a  number  of  elements.  Only  the 
location  of  comer  nodes  and  the  number  of  elements  along  the  PAFBLOCK  sides  need 
to  be  specified,  in  order  for  proportionately  sized  internal  elements  to  be  created. 
Typically  for  the  analysis  of  the  example  problems  described  in  Section  4,  there  were 
two  elements  per  PAFBLOCK  along  the  moving  boundary. 


3.2  Smoothing  of  nodal  movements 

In  the  course  of  preliminary  analyses  it  was  found  that  localised  effects  in  some  cases 
dismpted  the  solution  process,  particularly  for  highly  refined  meshes.  In  these  cases 
botmdary  node  stresses  were  found  to  be  highly  sensitive  to  the  positions  of 
neighboixring  boundary  nodes.  For  example,  a  node  protruding  slightly  from  the 
boundary  had  a  very  low  stress  and  was  hence  moved  too  far  outwards  (away  from 
the  hole  centre)  at  the  next  iteration.  The  neighbouring  nodes  would  then  be  left 
protmding  and  be  moved  even  further  in  the  next  iteration.  The  result  of  this  is  that 
alternate  nodes  would  Teap-frog'  each  other  by  increasing  amounts,  eventually 
causing  fatal  mesh  distortion.  To  circumvent  this  potential  problem  for  all  analyses, 
nodal  positions  as  defined  by  displacements  given  in  Equation  1,  were  adjusted  using 
three  point  averaging  to  include  the  effect  of  the  two  neighbouring  nodal  positions  as 
follows; 
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xf  =(x,._, +x,. +x,.^i)/3  (3) 

y-  =(y!-,+yi+yM)^^ 

where  x  ,  y  refer  to  the  nodal  co-ordinate  positions  determined  using  displacements 
given  by  Equation  (1),  and  x%y“  refer  to  the  values  after  smoothing. 


4.  Numerical  Examples  Relating  to  Holes  in  Plates 

In  this  section  the  results  obtained  from  application  of  the  method  to  sample  problems 
involving  stress  concentrators  in  two  dimensional  plates  are  presented.  For  all  analyses 
a  linear  elastic  constitutive  model  was  used  with  Yoxmg  s  modulus  E  -  70  000  MPa, 
and  Poisson's  ratio  v  =  0.32.  For  aU  cases  rework  optimisation  analyses  were 
undertaken  where  only  material  removal  occurs.  Fience  all  nodal  movements  were 
constrained  such  that  d.<0.  Also,  in  each  case  unless  otherwise  noted,  the  threshold 
value  of  stress  used  in  Equation  (1)  was  defined  as  ,  where 


i=k 


For  all  cases  involving  automated  botmdary  nodal  movement,  three  point  smoothing 
was  used  as  described  in  section  3.2. 


4.1  A  square  hole  in  a  2:1  biaxial  stress  field 

For  this  analysis  the  plate  geometry  as  shown  in  Figure  4  was  considered.  Here  the 
large  square  plate  contains  a  square  centrally  located  hole  with  20  mm  long  sides.  A 
remote  2:1  biaxial  stress  field  is  appHed  to  the  plate,  where  the  applied  stresses  are  50 
MPa  along  the  edges  x=±100  mm  and  100  MPa  along  the  edges  y=±100  mm.  Hence  one 
feature  of  this  initial  geometry  (with  the  sharp  corner)  is  that  significant  changes  in  the 
shape  are  expected,  and  therefore  is  considered  a  good  test  problem. 

Due  to  symmetry,  only  one  quarter  of  the  plate  was  modelled  in  the  finite  element 
analysis,  and  the  initial  mesh  discretisation  using  4  noded  elements  is  shown  in  Figure 
5.  The  optimisation  procedure  detailed  in  Section  3  was  used  with  the  parameters  r  = 
10  mm,  and  5  =  0.3.  The  solution  geometry  obtained  after  25  iterations  is  shown  in 
Figure  6.  It  can  be  seen  that  the  bormdary  is  positioned  entirely  outside  the  initial 
square  opening  as  required.  As  expected,  the  final  shape  is  the  known  solution  of  an 
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ellipse  with  an  aspect  ratio  of  two  to  one  for  the  major  and  minor  axes  [14].  For  this 
solution  the  tangential  stresses  along  the  boimdary  were  essentially  constant  with  a 
range  of  152.2  to  159.2  MPa,  noting  that  the  peak  stress  for  the  initial  geometry  was  535 
MPa  at  the  location  Ji:=y=10  mm.  The  contours  given  in  Figure  6  show  the  distribution 
of  maximum  absolute  principal  stress,  which  on  the  hole  boimdary  is  equivalent  to  the 
tangential  stress. 

To  determine  the  influence  of  the  value  of  the  step  factor  s ,  the  analyses  were  repeated 
taking  different  values  of  this  parameter.  The  relevant  results  are  given  in  Figure  7 
where  the  objective  function  versus  iteration  number  is  plotted  for  four  different  step 
factors.  It  can  be  seen  that  a  step  factor  of  0.4  causes  the  solution  to  fail  after  4  iterations 
when  the  boimdary  over  shoots  the  solution.  The  correct  shape  is  obtained  for  step 
factors  smaller  than  this  with  obviously  more  iterations  required  as  the  step  factor  is 
reduced.  At  a  step  factor  of  0.35  a  successful  solution  is  achieved  but  failure  almost 
occurs  at  iteration  5.  Robust  successful  solutions  are  obtained  at  values  of  0.3  and  0.2, 
where  it  is  seen  that  the  objective  function  decreases  monotonically  with  iteration 
number. 


4.2  A  round  hole  in  a  4:1  biaxial  stress  field 

For  this  analysis  the  plate  geometry  as  shown  in  Figure  8  was  considered.  Here  the 
large  rectangiilar  plate  contains  a  circular  centrally  located  hole  of  20  mm  diameter.  A 
remote  4  : 1  biaxial  stress  field  is  appHed  to  the  plate,  where  the  applied  stresses  are  50 
MPa  along  the  edges  x=±100  mm  and  200  MPa  along  the  edges  y=±100  mm.  One 
feature  of  the  initial  geometry  is  that  there  is  a  region  of  compressive  boimdary  stress, 
unlike  the  previous  example. 

Due  to  symmetry,  only  one  quarter  of  the  plate  was  modelled  in  the  finite  element 
analysis,  and  the  initial  mesh  discretisation  is  shown  in  Figure  9.  To  improve  accuracy 
a  highly  refined  mesh  consisting  of  8  noded  isoparametric  quadrilaterals  was  used. 
The  optimisation  procedure  detailed  in  Section  3  was  applied  with  the  parameters  r  = 
10  mm  and  .s  =0.1.  The  solution  geometry  obtained  after  250  iterations  is  shown  in 
Figure  10.  As  expected  the  final  shape  is  very  close  to  the  known  solution  of  an  ellipse 
with  an  aspect  ratio  of  four  to  one  for  the  major  and  minor  axes  [14],  which  is 
positioned  entirely  outside  the  circular  opening  as  required.  The  tangential  stress 
distribution  around  the  hole  boundary  is  shown  in  Figure  11.  For  this  solution  the 
tangential  stresses  along  the  boimdary  were  essentially  constant  (except  for  one  nodal 
position)  with  a  stress  range  of  250  to  275  MPa,  noting  that  the  peak  stress  for  the 
initial  geometry  was  550  MPa  at  the  location  x=10  mm  and  y=0  mm.  For  the  selected 
step  factor,  it  can  be  seen  that  the  optimisation  process  was  relatively  slow  to  converge 
(particularly  after  the  50th  iteration)  as  shown  in  Figure  12.  One  reason  for  this  is 
believed  to  be  due  to  the  fact  that  the  sensitivity  between  the  length  of  the  ellipse  and 
the  boimdary  stress  is  very  low.  Along  the  sides  of  the  eUipse,  very  small  nodal 
movements  are  required  to  find  the  optimal  solution  and  there  is  a  high  level  of 
dependence  between  the  nodes.  For  the  last  100  iterations  it  was  necessary  to  turn  off 
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the  three  point  smoothing  as  it  was  prone  to  stifle  the  very  small  nodal  movements. 
From  the  monotonically  reducing  objective  function  shown  in  Figure  12,  it  is  clear  that 
use  of  a  larger  step  factor  would  have  reduced  the  number  of  iterations  required.  The 
compressive  stress  region  in  the  mitial  configuration  did  not  effect  the  performance. 


4.3  Corrosion  removal  in  a  uniaxial  stress  field 

The  size  of  full  depth  skin  corrosion  removal  in  aircraft  structures  is  often  restricted  by 
neighbouring  structural  details.  In  the  present  analysis  the  optimisation  procedure  is 
applied  to  a  problem  which  addresses  this  issue.  For  this  analysis  the  plate  geometry 
as  shown  in  Figure  13  was  considered.  Here  the  large  rectangular  plate  contains  a 
centraUy  located  circular  hole  of  12  mm  diameter.  A  remote  uniaxial  stress  field  of  100 
MPa  is  applied  to  the  plate  along  the  edges  y=±100  mm.  The  key  feature  is  that 
allowable  material  removal  is  constrained  by  the  lines  at  y=±12  mm. 

Due  to  synunetry,  only  one  quarter  of  the  plate  was  modelled  in  the  fmite  element 
analysis,  and  the  initial  mesh  discretisation  is  shown  in  Figure  14.  Separate  analyses 
were  conducted  using  either  all  four  noded  or  all  eight  noded  elements.  A  manual 
implementation  of  the  moving  boimdary  method  described  in  Section  3  was  used.  The 
optimisation  was  performed  by  manually  moving  PAFBLOCK  comer  nodes  based  on 
a  chosen  threshold  stress  and  the  values  of  boundary  tangential  stress.  Where 
boundary  hoop  stress  was  below  the  threshold  value,  nodes  were  moved  inwards  to 
the  threshold  stress  contour  of  maximum  principle  stress.  The  new  coordinates  were 
edited  into  the  PAFEC  data  file  of  the  previous  iteration  The  solution  geometry 
obtained  after  14  iterations  is  shown  in  Figure  15.  It  can  be  seen  that  a  uniform 
tangential  stress  occurs  over  the  curved  region  of  the  boxmdary,  which  is  also  shown 
more  clearly  in  Figure  16.  The  peak  stress  on  the  hole  boundary  was  initially  311  MPa 
which  has  been  reduced  to  essentially  a  constant  value  along  the  curved  part  of  the 
boundary,  in  the  range  167  to  178  MPa.  There  is  good  agreement  between  the  analyses 
using  either  all  4  noded  or  all  8  noded  elements,  although  as  expected  the  8  noded 
solution  is  better  (i.e.  less  deviation  from  the  mean  along  the  curved  boxmdary). 

This  example  highlights  the  disadvantages  of  using  Von  Mises  stress  as  a  criterion  for 
removing  material,  as  used  by  some  authors.  Because  of  its  root-mean-square 
formulation,  boxmdary  regions  with  compressive  tangential  stresses  will  have  positive 
Von  Mises  values.  Hence  for  such  situations  the  adoption  of  the  Von  Mises  based 
removal  criteria  would  erroneously  not  remove  material  in  these  regions. 

At  the  level  of  mesh  refinement  used,  the  abmpt  comer  shown  in  Figure  15  is  free  of 
raised  stresses  because  of  the  unidirectional  remote  loading  and  the  optimal  rework 
shape.  Any  transverse  load  component  would  have  led  to  a  roxmded  shape  in  this 
comer  region.  It  is  interesting  to  note  that  the  curved  part  of  the  boxmdary  is  neither 
elliptical  or  circular  in  form  and  is  xmhkely  to  be  obtained  by  analytical  methods. 
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4.4  Inclined  oval  hole  in  a  uniaxial  stress  field  (representative  of  F-111 
fuel  flow  vent  hole  number  13) 

As  discussed  in  the  previous  example  the  shape  of  a  rework  is  often  constrained  by 
geometric  limits.  This  issue  is  further  explored  in  the  present  example.  For  this  analysis 
the  plate  geometry  as  shown  in  Figure  17  was  considered.  Flere  the  large  rectangular 
plate  contains  a  centrally  located  oval  shaped  hole  inclined  at  16  de^ees  to  the 
direction  of  uniaxial  loading.  The  remote  uniaxial  stress  field  of  100  MPa  is  applied  to 
the  plate  along  the  edges  y=±140  mm.  The  key  feature  of  the  analysis  is  that  material 
removal  is  constrained  by  the  line  shown  in  Figure  17,  which  represents  a  severe 
geometric  constraint.  This  geometry  has  been  chosen  to  be  approximately 
representative  of  the  fuel  flow  vent  hole  nxunber  13  (FFVH13)  region  in  the  stiffener 
near  the  upper  wing  skin  of  the  F-111  wing  pivot  fitting  (which  is  a  well  known  crack 
initiation  site  [15]).  The  plate  geometry  shown  in  Figure  17  is  a  two  dimensional 
idealisation  representing  the  stiffener,  while  the  location  of  the  constraint  line 
corresponds  to  the  location  of  the  upper  wing  skin,  with  its  thickness  in  the  out  of 
plane  direction.  In  recent  years  the  FFVH13  region  has  been  subjected  to  grind-out 
rework  profiles  to  remove  service  cracks  and  reduce  the  extent  of  the  stress 
concentration  to  some  extent  [15]. 

The  initial  mesh  discretisation  for  this  problem,  consisting  of  four  noded  elements,  is 
shown  in  Figure  18.  The  optimisation  procediue  detailed  in  Section  3  was  applied  with 
the  foUowing  parameters:  (i)  ex,,,  =  0.5(o-„,  +0’^’')/  (h)  r  =  6.35  mm,  and  (iii)  5=  0.3. 
Iterations  were  continued  until  the  reduction  in  tangential  stress  at  point  A  became 
less  than  1%.  The  geometry  at  the  final  iteration  is  shown  in  Figure  19,  while  the 
associated  distribution  of  tangential  stress  around  the  hole  boxmdary  is  shown  in 
Figure  20.  It  can  be  seen  that  the  process  provided  a  good  result  in  the  critical  region  of 
point  A.  Here  the  initial  peak  stress  of  309  MPa  was  reduced  to  a  value  of  171  MPa 
which  was  rendered  uniform  in  the  local  region  (i.e.  nodes  11  to  27).  However  near 
point  B  (nodes  37  and  42)  it  can  be  seen  that  the  presence  of  the  constraint  line  has 
reduced  the  benefit  of  the  optimisation  process.  Away  from  the  constraint  line  the  final 
shape  is  approximately  elliptical  in  form.  This  example  has  demonstrated  clearly  that 
shape  optimisation  results  are  highly  dependent  on  the  direction  of  the  surrounding 
principle  stresses,  in  conjunction  with  the  location  of  geometric  constraints.  It  has  also 
shown  that  near  point  B  the  initial  geometry  was  close  to  optimal. 


5.  Numerical  Examples  Relating  to  Bonded 
Reinforcements 


Bonded  reinforcements  and  repairs  to  engineering  components  function  primarily  by 
transferring  load  from  the  component  to  the  patch  by  shear  deformation  of  the 
adhesive.  This  reduces  the  magnitude  of  strains  in  the  repaired  component  and 
typically  results  in  increased  fatigue  hves  and  static  strength.  One  issue  concerning  the 
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improved  design  of  bonded  repairs  which  needs  to  addressed  is  the  existence  of  large 
peaks  in  the  adhesive  shear  strain  distribution.  For  example,  a  severe  peak  typically 
occurs  near  the  end  of  the  patch.  Such  peaks  have  the  potential  to  cause  failure  of  the 
adhesive  system  and  compromise  the  performance  of  the  repair.  In  this  section  the 
method  explained  in  section  3  is  slightly  modified  to  extend  its  application  to  the 
optimisation  of  the  shape  of  the  bonded  reinforcement,  while  minimising  the  adhesive 
stress  levels.  Two  example  problems  were  considered,  where  in  each  case  a  two 
dimensional  analysis  was  undertaken  using  a  linear  elastic  constitutive  model.  In  the 
first  case  the  patch  thickness  was  varied  continuously  while  in  the  second,  for  a  given 
patch  design,  the  adhesive  thickness  was  varied  continuously. 


5.1  Bonded  joint  taper 

The  configuration  under  study  was  a  double  lap  joint  as  shown  in  Figure  21.  It  should 
be  noted  that  this  geometry  is  representative  of  a  two  dimensional  idealisation  of  a 
bonded  repair  to  a  cracked  plate.  The  inner  adherend  is  a  plate  of  4  mm  thickness 
which  is  loaded  by  a  remote  stress  of  100  MPa.  An  outer  adherend  (i.e.  patch)  is 
bonded  to  each  side  of  this  plate  by  an  adhesive  layer  having  a  uniform  thic^ess  of 
0.15  mm.  Thereby,  symmetry  is  retained  with  respect  to  the  plate  mid-plane  {y=0),  and 
hence  plate  bending  is  eliminated.  The  thickness  of  each  patch  is  denoted  to  and  is 
2  mm.  As  is  commonly  advocated,  an  initial  linear  taper  with  a  1:10  slope  was  used  at 
the  ends  of  each  patch  to  reduce  the  magnitude  of  the  adhesive  stress  concentration  at 
the  end  of  the  patch.  For  this  analysis  a  linear  elastic  constitutive  model  was  used  to 
model  the  plate  and  the  patch  with  Young's  modulus  E  =  70  000  MPa,  and  Poisson's 
ratio  V  =  0.32  In  this  optimisation  we  wish  to  change  the  shape  of  the  tapered  region  to 
achieve  a  constant  adhesive  stress  distribution.  The  procedure  used  was  to  reduce  the 
patch  thickness  (at  a  given  x  co-ordinate)  where  the  adhesive  stresses  were  high,  while 
increasing  patch  thickness  where  adhesive  stresses  were  low.  Hence  a  slight  variation 
of  Equation  (1)  to  represent  movement  of  nodes  on  the  free  boundary  of  the  patch  was 
used  as  follows 


=  (6) 

'^th 


where  positive  d.  indicates  material  addition  from  the  patch  outer  botmdary,  r  -  is  the 
shear  stress  at  the  adhesive  midplane  node  i ,  is  a  threshold  shear  stress  at  the 
adhesive  midplane,  and  r  and  s  are  as  previously  defined. 

Due  to  symmetry,  only  one  quarter  of  the  joint  was  modelled  in  the  fmite  element 
analysis,  and  the  initial  mesh  discretisation  is  shown  in  Figure  22(a).  The  mesh 
consisted  of  4  noded  quadrilateral  elements.  A  very  dense  mesh  was  used  near  the  end 
of  the  patch  to  accurately  model  the  high  stress  gradients  in  this  region.  For  the 
adhesive  a  Young's  Modulus  of  840  MPa  and  a  Poisson  s  Ratio  of  0.3  was  used.  Other 
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material  properties  were  as  per  section  5.1.  The  optimisation  procedure  was  applied 
with  the  parameters  r  =  2  mm  and  5=0.1.  One  further  practical  constraint  was 
introduced  such  that  the  patch  thickness  could  not  be  less  that  0.01  mm  and  not  greater 
than  2  mm.  The  solution  geometry  obtained  after  89  iterations  is  shown  in  Figure  22(b). 
As  a  comparison  a  one  dimensional  theoretical  analysis  given  in  [16]  was  formd  to  give 
the  profile  shown  in  Figure  22(c).  The  equation  from  [16]  defining  the  thickness  in  the 
tapered  region  as  a  function  of  position  is 


where  xis  the  distance  from  the  start  of  the  tapered  region,  is  the  maximum 
thickness  of  the  patch,  and  I  is  the  length  of  the  tapered  region. 

In  Figure  23  the  adhesive  shear  stresses  are  plotted  for  the  following  three  taper 
profiles:  (i)  linear  taper,  (ii)  nrrmerically  optimised  taper,  and  (iii)  analytically 
optimised  taper  [16].  It  can  be  seen  that  the  process  provided  a  good  result  in  reducing 
the  peak  shear  stress,  and  rendering  the  shear  stresses  relatively  uniform.  Here  the 
initial  peak  stress  of  8.2  MPa  for  the  linear  taper  case  has  been  reduced  by  30%.  The 
one  dimensional  analysis  gives  a  similarly  significant  reduction  in  peak  adhesive 
stress.  It  is  included  here  as  a  convenient  comparison  for  this  particular  loading  case.  It 
must  be  noted  that  the  advantage  of  the  finite  element  method  is  that  key  parameters 
such  as  loading  t5q)e  and  material  properties  are  easily  changed,  and  their  modification 
is  expected  to  pose  no  particular  difficulty. 

It  is  interesting  to  note  that  for  the  given  inner  adherend,  a  reduction  in  the  adhesive 
peak  at  the  end  the  patch  can  only  be  achieved  through  changing  one  or  more  of  the 
material  or  geometric  properties  of  the  adhesive  or  the  patch,  near  the  end  of  the  patch 
[17].  For  example  by  increasing  the  adhesive  thickness.  While  the  optimisation  process 
has  been  demonstrated  for  the  case  of  rninirnising  adhesive  shear  stress,  it  is  believed 
other  stress  quantities,  such  as  von  Mises  stress  could  have  been  selected. 


5.2  Bonded  repair  to  F/A-18  470  bulkhead 

For  this  problem  a  semi-symmetric  two  dimensional  MSC/NASTRAN  finite  element 
model  obtained  from  reference  [18]  was  used  as  a  starting  point.  The  relevent  finite 
element  geometry  and  mesh  is  shown  in  Figure  24.  The  location  at  Point  A  is  a  known 
crack  initiation  site  for  the  bulkhead  rmder  the  action  of  flight  loads,  due  to  the  stress 
magnification  in  this  local  region.  One  life  extension  option  under  consideration  at 
AMRL  for  reducing  the  magnitude  of  the  stress  concentration  involves  bonding  an  8 
ply  boron  patch  to  the  critical  region  of  the  bulkhead.  For  convenience  a  PAFEC  finite 
element  model  was  developed  to  represent  the  substructured  critical  local  region  as 
shown  in  Figure  25.  Here  it  can  be  seen  that  extra  material  was  added  to  the  right  of 
the  substructured  region  to  enable  correct  distribution  of  local  stresses.  Hence  the 
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magnitudes  of  the  two  uniaxial  loads  (342386  N  and  42675  N)  were  chosen  such  that 
the  stresses  along  the  line  of  symmetry  were  the  same  as  those  given  in  the  full 
MSC/NASTRAN  model. 

The  finite  element  mesh  for  the  substructured  region  with  the  8  ply  boron/  epoxy 
patch  is  shown  in  Figure  26.  For  this  analysis  the  boron  patch  was  treated  as  an 
isotropic  material  with  a  Young's  Modulus  of  213,000  MPa  and  a  Poisson's  Ratio  of  0.3. 
The  material  properties  of  the  plate  and  adhesive  were  as  used  in  Section  5.1.  Here  the 
adhesive  thickness  was  constant  at  0.13  mm  and  the  patch  thickness  was  1.04  mm.  The 
Von  Mises  stress  results  for  the  adhesive  are  given  in  Figure  29,  with  the  associated 
reduction  in  stresses  at  point  A  in  the  bulkhead  being  34%  as  compared  to  the 
unpatched  results.  From  Figure  29  it  can  be  seen  that  the  adhesive  Von  Mises  stresses 
vary  along  the  joint  and  peak  at  point  B.  Clearly,  an  improved  repair  design  would 
require  a  more  uniform  adhesive  stress  distribution,  with  an  associated  reduced  peak 
value. 

To  optimise  the  adhesive  stress  distribution,  the  same  procedure  as  given  in  Section  5.1 
was  used,  except  in  this  case  the  thickness  of  the  adhesive  layer  was  varied.  This  was 
achieved  by  moving  nodes  on  the  interface  between  the  adhesive  and  the  patch  using 
the  relation  given  in  equation  (6)  with  slight  modification.  Here  the  shear  stress  terms 
are  replaced  by  Von  Mises  stress,  and  the  adhesive  thickness  is  increased  at  high  stress 
regions  and  reduced  at  low  stress  regions.  For  this  analysis  the  location  of  the  outer 
boundary  of  the  patch  remains  fixed  for  all  iterations.  The  initial  geometry  of  the 
adhesive  layer  consisted  of  a  xmiform  thickness  of  0.13  mm.  While  for  the  patch,  the 
initial  geometry  was  selected  to  have  the  following  outer  profile  (i)  nodes  1  to  5  were 
on  the  arc  of  a  circle  of  radius  25  mm,  and  nodes  5  to  12  were  on  a  straight  line.  Hence 
the  patch  had  a  variation  in  cross  sectional  thickness  in  the  y  direction,  with  the 
thicknesses  at  nodes  1, 3, 5  and  12  being  3.12, 2.92, 2.67  and  1.39  mm  respectively.  The 
optimisation  procedure  was  applied  with  the  parameters  r  =  1.5  mm  and  s  —0.1.  One 
further  practical  constraint  was  introduced  such  that  the  adhesive  thickness  could  not 
be  less  tiian  0.13  mm. 

The  bonded  repair  geometry  shown  in  Figure  27  was  achieved  after  47  iterations.  This 
profile  has  a  much  greater  stress  reduction  at  point  A  (42%)  and  a  more  uniform 
adhesive  stress  distribution  with  a  lower  peak  value,  as  compared  to  the  standard 
patch  design.  To  enable  a  better  comparison  with  the  nominal  8  layer  patch  case  the 
iteration  process  was  continued  imtil  the  same  stress  reduction  of  34%  (after  a  further 
53  iterations)  was  achieved  at  point  A  (in  this  process  the  adhesive  layer  becomes 
thicker).  The  final  bonded  repair  geometry  shown  for  this  optimisation  analysis  is 
shown  in  Figure  28,  and  the  adhesive  stresses  are  given  in  Figure  29.  It  can  be  seen  that 
the  adhesive  stresses  have  been  further  reduced  and  remain  uniform.  As  can  be 
expected,  these  results  are  consistent  with  the  trend  that  generally  in  bonded 
reinforcement  analyses,  a  lower  patch  effectiveness  is  associated  with  a  lower  adhesive 
stress. 
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6.  Discussion 


6.1  Typical  results 

The  example  problems  that  have  been  presented  (along  with  those  in  the  prominent 
hterature)  show  that  there  are  often  large  benefits  to  determining  the  optimal  shape  for 
a  given  structural  detail  imder  specified  loading  conditions.  Reductions  in  peak  stress 
of  30%  and  upwards  can  be  achieved  by  accurate  shaping  of  the  boimdary.  The  use  of 
circular  form  shapes  in  fatigue  prone  structures  is  a  poor  alternative  to  the  optimal  free 
form  shape  in  cases  where  the  principal  stress  directions  are  known  not  to  vary  with 
load  case.  For  many  problems  the  optimal  solution  is  characterised  by  tmiform 
boxmdary  stresses  in  the  critical  region,  as  discussed  in  Section  2.2  and  demonstrated 
by  the  results  presented  in  Sections  4  and  5. 

In  relation  to  optimisation  of  bonded  repairs  the  problem  of  optimal  shaping  away 
from  a  cracked  region  is  solved.  Examples  are  presented  with  good  results  showing 
two  distinct  options.  In  the  first  option  the  patch  thickness  is  varied  continuously  to 
provide  a  reduced  and  uniform  adhesive  stress  distribution  while  still  maintaining  the 
same  load  transfer  to  the  patch.  In  the  second,  for  a  given  patch  design,  the  adhesive 
thickness  is  varied  continuously  to  give  a  reduced  and  uniform  adhesive  stress 
distribution  while  still  effectively  reducing  stress  in  the  reinforced  component. 


6.2  Numerical  implementation  issues 

It  is  beheved  that  the  method  developed  here  is  a  simple  approach  that  can  be  adapted 
to  suit  finite  element  systems  that  do  not  provide  an  optimisation  feature.  For  example 
the  software  developed  to  interface  with  the  PAFEC  finite  element  code  (Appendix  A) 
contams  only  several  hundred  lines  of  code  with  minimal  branching. 

In  some  of  the  example  problems  highly  refined  meshes  were  used  which  clearly 
increased  the  number  of  iterations  required.  This  was  not  considered  a  concern  since 
the  total  solution  time  was  fast.  However  if  solution  time  is  a  concern,  one  obvious 
tactic  is  not  to  overly  refine  the  mesh.  Furthermore,  consideration  can  also  be  given  to 
using  a  coarse  mesh  initially  and  then  changing  to  a  more  refined  version  when  the 
optimal  solution  is  obtained.  Similarly  no  concerted  attempt  was  made  to  minimise  the 
number  of  iterations  by  changing  the  parameters  in  Equation  (1),  again  because  of  the 
fast  solution  time.  Nevertheless  if  desired,  it  is  beheved  that  there  is  significant  scope 
to  change  these  key  parameters  without  degrading  the  accuracy  of  the  solution. 
Similarly  while  a  simple  three  point  smoothing  method  was  successfully  apphed  to 
determine  nodal  displacements,  a  more  sophisticated  smoothing  method  could  readily 
be  used  if  desired. 

The  appropriate  choice  of  the  threshold  stress  in  Equation  (1)  is  considered  a  powerful 
feature  allowing  either  material  removal  or  addition,  or  a  weighted  combination  of 
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both.  It  is  also  interesting  to  note  that  because  of  the  formulation  of  Equation  (1),  the 
absolute  magnitude  of  bormdary  movement  steps  decreases  automatically  with 
increasing  iteration  number,  as  the  optimal  solution  is  approached.  Hence  the  use  of  a 
uniform  step  factor  was  foimd  to  be  very  effective. 


7.  Conclusion 

This  report  presents  the  development  and  automated  numerical  implementation  of  a 
gradientless  optimisation  method  for  the  analysis  of  problems  relating  to  life  extension 
of  aircraft  components,  namely  (i)  shape  reworking  of  stress  concentrators  and  (ii) 
shape  changes  for  bonded  reinforcements.  The  method  has  been  implemented  to  be 
used  in  conjunction  with  the  finite  element  code  PAFEC,  and  it  is  considered  suitable 
for  use  with  other  typical  standard  commercial  finite  element  packages.  The  key 
feature  of  the  approach  is  to  achieve  constant  boundary  stresses  in  regions  of  interest, 
by  correctly  moving  boundary  nodal  positions.  Typically  nodes  on  the  stress 
concentrator  boundary  are  moved  by  an  amount  dependent  on  the  sign  and 
magnitude  of  the  local  hoop  stress  (relative  to  a  selected  threshold  value)  obtained 
from  a  previous  iteration  of  a  standard  finite  element  analysis.  The  results  of  the 
numerical  examples  presented  demonstrate  the  effectiveness  of  the  approach  that  has 
been  developed. 

For  the  example  rework  problems  that  have  been  presented,  along  with  those  in  the 
prominent  hterature,  it  has  been  shown  that  there  are  often  large  benefits  brought 
about  by  finding  the  optimal  shape  for  a  given  structural  detail.  Reductions  in  peak 
stress  of  30%  and  upwards  can  be  achieved  by  accurate  shaping  of  the  boundary 
geometry.  The  use  of  circular  or  even  elliptical  form  shapes  in  fatigue  prone  structures 
is  a  poor  alternative  to  the  optimal  free  form  shape.  In  relation  to  optimisation  of 
bonded  repairs  the  problem  of  optimal  shaping  away  from  a  cracked  region  has  been 
solved.  Examples  are  presented  with  good  results  showing  two  distinct  options.  In  the 
first  option  the  patch  thickness  is  varied  continuously  to  provide  a  reduced  and 
uniform  adhesive  stress  distribution  while  still  maintaining  the  same  load  transfer  to 
the  patch.  In  the  second,  for  a  given  patch  design,  tiie  adhesive  thickness  is  varied 
continuously  to  give  a  reduced  and  uniform  adhesive  stress  distribution  while  still 
effectively  reducing  stress  in  the  reinforced  component. 

It  is  considered  that  a  number  of  extensions  to  the  work  are  warranted.  Firstly,  it  is 
believed  that  there  is  significant  scope  to  change  key  parameters  in  the  algorithm 
procedure  to  reduce  the  number  of  iterations  required,  if  desired.  Secondly,  the 
approach  seems  directly  appHcable  to  three  dimensional  geometries. 
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Figure  3.  Flowchart  for  implementation  of  moving  boundary  method  with  standard  finite 
element  analysis. 
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Figure  4.  Geometry  and  loading  arrangement  for  a  large  square  plate  containing  a  square 
hole  with  a  remote  2:1  biaxial  stress  field. 


Figure  5.  Finite  element  mesh  for  modelling  a  large  square  plate  containing  a  square  hole. 
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Figure  6.  Final  shape  far  problem  of  optimisation  of  an  initial  square  hole  in  a  large  plate 

vnth  a  remote  2:1  biaxial  stress  field,  showing  distribution  of  maximum  principal  stress. 
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Figure  7. 


Convergence  of  objective  function  for  various  step  factors  for  a  large  plate 
containing  an  initial  square  hole  with  a  remote  2:1  biaxial  stress  field. 
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Figure  8.  Geometry  and  loading  arrangement  for  a  large  square  plate  containing  a  round 
hole  with  a  remote  4:1  biaxial  stress  field. 


Figure  9.  Finite  element  mesh  for  modelling  a  large  square  plate  containing  a  round  hole. 
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Figure  10.  Final  hole  shape  for  problem  of  optimisation  of  an  initial  round  hole  in  a  large  plate 
with  a  remote  4:1  biaxial  stress  field. 
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Figure  11.  Distribution  of  tangential  stress  around  hole  boundary  of  final  shape  for  initial 
round  hole  in  a  large  plate  with  a  remote  4:1  biaxial  stress  field 


Figure  12.  Convergence  of  objective  function  for  a  large  plate  containing  an  initial  round  hole 
with  a  remote  4:1  biaxial  stress  field. 
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Figure  14.  Finite  element  mesh  for  modelling  a  large  square  plate  containing  a  round  hole. 
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Figure  15.  Final  shape  for  problem  of  constrained  optimisation  of  an  initial  round  hole  in  a 

large  uniaxially  loaded  plate,  showing  distribution  of  maximum  principal  stresses. 


Figure  16.  Distribution  of  tangential  stress  around  hole  boundary  of  final  shape  for  initial 
round  hole  in  a  large  uniaxially  loaded  plate. 
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Figure  19.  Final  shape  for  problem  of  constrained  optimisation  of  an  initial  inclined  oval  hole 
in  a  large  uniaxially  loaded  plate,  showing  distribution  of  maximum  principal 
stresses. 


Figure  20.  Distribution  of  tangential  stress  around  hole  boundary  of  final  shape  for  constrained 
optimisation  of  an  initial  inclined  oval  hole  in  a  large  uniaxailly  loaded  plate. 
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Figure  22(a).  Finite  element  mesh  for  bonded  double  lap-joint  with  initial  linear  taper 
profile  for  patch  (only  region  near  taper  shown). 


Figure  22(b).  Final  shape  for  optimisation  of  the  tapered  region  of  patch  for  bonded 
double  lap-joint  (only  region  near  taper  shown). 
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Figure  23.  Adhesive  shear  stress  distributions  for  a  bonded  double  lap-joint  with  a 

unidirectional  applied  stress,  for  various  taper  profiles  (in  taper  region  only). 
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Figure  26.  Finite  element  mesh  for  critical  region  ofF/A-18  470  bulkhead  substructure  for 
nominal  8  layer  bonded  patch  with  stress  reduction  at  point  A  of  34% . 


Figure  27.  Finite  element  mesh  for  critical  region  ofF/A-18  470  bulkhead  substructure  for 
optimised  patch  with  stress  reduction  at  point  A  of  42%,.). 
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Figure  28.  Finite  element  mesh  for  critical  region  ofF/A-18  470  bulkhead  substructure  for 
optimised  patch  with  stress  reduction  at  point  A  of  34%. 


Nodes  around  boundary 
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Figure  29.  Comparison  of  von  Mises  adhesive  stress  distributions  for  F/A-18  470  bulkhead 
substructure  region  with  initial  and  optimised  bonded  patches. 
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9.  Appendix  A: 

Implementation  of  Moving  Boundary  Method  In 
PAFEC  Finite  Element  Code 


9.1  General  comments 

The  implementation  of  the  moving  boundary  node  method  with  PAFEC  is  described 
by  the  flow  chart  given  in  Figure  A.l.  The  finite  element  software  is  used  to  solve  for 
the  stresses  at  each  iteration  with  the  finite  element  model  constructed  such  that  the 
boundary  nodes  to  be  moved  are  comers  of  PAFEC  PAFBLOCKS.  PAFBLOCKS  are  a 
geometric  entity  that  contain  proportional  iso-meshing  of  finite  elements  within.  A 
new  PAFBLOCK  comer  position  leads  to  new  positions  of  all  the  internal  nodes  that 
maintain  the  original  mesh  proportions.  It  is  advantageous  in  some  cases  to  select  the 
initial  shapes  of  the  PAFBLOCKS  such  that  the  new  comer  positions  will  not  cause  any 
elements  to  become  badly  distorted  as  the  iteration  process  proceeds. 

The  initial  model  is  prepared  in  the  normal  way  using  the  operating  system  and  a  text 
editor.  Overall  aspects  of  the  process  such  as  the  iterations  and  file  handling  are 
controlled  by  a  UNIX  sheU  script  which  is  given  in  section  9.2.  Convergence  checking 
and  node  movements  are  performed  by  a  FORTRAN  program  which  is  executed  by 
the  shell  script.  The  FORTRAN  program  needs  to  be  modified  sHghtly  for  different 
problems.  Changes  mostly  relate  to  the  method  used  to  calculate  the  direction  of  node 
movements,  and  three  point  smoothing  for  various  levels  of  symmetry. 

A  dedicated  directory  is  preferred  for  each  problem  due  to  the  amount  of  file  handling 
performed  by  the  shell  script  and  the  PAFEC  solution.  The  primary  output  is  a  log  file 
containing  mainly  stresses  and  nodal  coordinates  for  each  iteration.  The  PAFEC  binary 
file  containing  the  mesh  and  model  description  is  also  retained  so  that  any  of  the  shape 
variations  can  be  plotted.  A  second  shell  script  is  used  to  delete  and  initialise  files  for 
restarting  the  problem. 
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Figure  A.l  -  Flow  chart  for  implementation  of  the  moving  boundary  node  method  with  PAFEC 

finite  element  code 
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9.2  Unix  shell  script  for  controlling  iterations 


maxcount="31 " 

read  j  obname  <  JNAME . DAT 
count="l" 

# 

until  [  $count  =  $maxcount  ] 
do 

echo  $count 
pa f del  $j obname 
pafrun  $j obname 

cp  $jobname" -BS"  $ jobname$count" .BS" 

pupname="PUP. JNL" 

echo  $pupname 

pup  $pupname 

mn 

# 

#  Terminate  job  on  Solution  solved 

# 

if  [  -f  JOB. SOLVED  ];  then 
break 
f  i 

count='expr  $count  +  1' 
done 


9.3  Fortran  program  for  moving  nodes 

program  movenodes 

C  This  is  the  version  of  code  used  for  the  analysis  presented  in  section  4.2. 
C  Here  it  is  assumed  that  the  problem  domain  is  in  the  first  quadrant  and  the 
C  plate  geometry  is  quarter  symmetric. 

C  Declarations 

implicit  none 
logical  solved 

character  aline*80, jobname*10, filename*20 
integer  eofn, i, j , k, ios, n, nodes, joblen, linelen 
real  rnum, x, y, z, si, s2, s3 , thl, th2, th3,pi 

real  hoop, hooptot, avhoop, r, th, d, h, dhoop, vartot, 

variance 

character  chr 

dimension  n(40,20) ,x{40,  20) ,y(40,20) , z (40,20) , hoop (40,20)  , 
th (40, 20) ,d (40,20) , si (40,20) ,s2 (40,20) ,s3 (40,20) , thl (40,20) , 

*th2 (40, 20) , th3 (40,20) , variance (20) , avhoop (20) 


C  Read  in  iteration  number 

open  (unit=10,file=' ITNUM.DAT') 
read(10,235)  j 
235  format (i2) 
close  (10) 


C  Open  PUPPIES  output  file 

f ilename=  'JOB . PUPPIES ' 

open  (unit=20, file=filename, err=1010, status=' old' ) 

C  Read  through  alphajunk  to  find  first  line  of  numeric  data 
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C  (i.e.  nodal  coordinates) 
ios=l 

do  while  (ios.ne.O) 
read (20, *, iostat=ios)  rnum 
end  do 

baclcspace  (20) 

C  Read  in  node  numbers  and  coordinates 
i=0 

do  while  (ios.eq.O) 

i=i+l  ^  ,  ,  , 

read (20, *, iostat=ios)  n(i,j),  x(i,j),  y{i,j)/  z(i,3) 
end  do 
nodes=i-l 


C  Read  through  more  alphajunk  to  find  first  line  of  stress  data 
ios=l 

do  while  (ios.ne.O) 
read (20, *, iostat=ios)  rnum 
end  do 

backspace (20) 

C  Read  in  stresses  and  directions 
i=0 

do  while  (ios.eq.O) 
i=i+l 

read (20, * , iostat=ios )  n(i,j),sl(i,j),s2(i,j), 

*  s3 (i, j) ,thl (i, j) , th2 (i, j) ,th3 (i, j ) 

end  do 
close (20) 

C  If  this  is  not  the  first  iteration  get  the  results  of  previous 
C  iterations  from  the  log  file  JOBOPT.LOG 
if  (j.gt.l)  then 
k=j 

open  (unit=10 , f ile= ' JOBOPT . LOG ' , status= ' old ' ) 
do  j=l,k-l 

read (10, *, iostat=ios)  rnum 
read  (10, *, iostat=ios)  avhoop(j) 
read  (10, *, iostat=ios)  variance ( j ) 
read (10, *, iostat=ios)  rnum 

do  i=l, nodes  _  _  ^  ,  ,  . 

read (10, *, iostat=ios)  n(i,j),x(i,j),y(i,j),sl(i,j),s2(i,3), 

*  thl(i, j) ,th2(i, j) ,th(i, j) ,d(i, j) 
end  do 

end  do 
j=k 

end  if 
close  (10) 

C  Store  hoop  stresses  in  array  hoop  and  calculate  average  hoop  stress 
hooptot=0 
do  i=l, nodes 

if  (abs (si (i, j ) ) .gt.abs (s2  (i, j ) ) )  then 
hoop  (i, j ) =sl  (i, j ) 
else 

hoop  (i, j ) =s2  (i, j ) 
end  if 

hooptot=hooptot+hoop (i,  j ) 
end  do 

avhoop  ( j ) =hooptot/nodes 

C  Calculate  variance 
vartot*0 
do  i=l, nodes 
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vartot=vartot+ (hoop(i, j ) -avhoop( j ) ) **2 
end  do 

variance (j ) =sqrt (vartot/avhoop ( j ) ) 
Calculate  base  length 

r=sqrt (x(l,j) **2+y(l,j)**2) 
r=6 


C  Calculate  the  direction  to  move  nodes  (clockwise  from  y  axis) 
pi=3. 14159 
do  i=l, nodes 
if  (i.eq.l)  th(i,j)=0 
if  (i. gt.l. and. i. It. nodes)  then 

h=sqrt( (x(i+l, j)-x(i-l, j) )**2+(y(i+l, j)-y(i-l,  j) )**2) 

th(i, j)=asin( (y (i+1, j ) -y (i-1, j ) ) /h) 
th(i, j)=-th(i, j) 
end  if 

if  (i.eq. nodes)  th(i,j)=pi/2 
end  do 

C  Now  calculate  d,  the  distance  to  be  moved  and  check  for  convergence 
solved=.true. 
do  i=l, nodes 

dhoop= (avhoop ( j ) -hoop (i, j ) ) /avhoop ( j ) 
d (i, j ) =dhoop*r*0 . 7 

if  (abs  (dhoop) .gt.0.02)  solved=. false, 
end  do 

C  Write  LOG  file  to  JOBOPT.LOG 
f i lename= ' JOBOPT . LOG ’ 
open  (unit=30, file=filename) 

k=j 

do  j=l,k 

write  (30,*)  j,'  =  iteration  number' 
write(30,*)  avhoop(j),'  =  average  hoop  stress' 
write(30,*)  variance (j ), '  =  Variance  ' 
write(30,*)  'NODE  X  Y  SI  S2'// 

*  '  thl  th2  th(i)  d(i)  ' 

do  i=l, nodes 

write (30,220)  n(i,j),x(i,j),y(i,j),sl(i,j),s2(i,j), 

*  thl (i, j) ,th2 (i,j),th(i,j),d(i,j) 

220  format (i3, 8f8 . 2) 

end  do 
end  do 
j=k 

close  (30) 

if  (solved)  then 
open (unit=30, file- 'JOB. SOLVED ' ) 
close (30) 
stop 
end  if 

C  Change  coordinates 
do  i=l, nodes 

if  (d (i, j ) .gt . 0)  then 
x(i, j)=x(i, j)+d(i, j)*sin(th(i, j) ) 
y(i, j)=y(i/ j)+d(i, j) *cos (th(i, j) ) 
end  if 
end  do 

C  Make  copy  of  pafec  data  file  in  TEMP.DAT 
filename- 'JOB. DAT' 

open  (unit-30,  file-filename, err-1 020, status-' old' ) 
open  (unit-40,  file- 'TEMP. DAT' ) 


37 


DSTO-RR-0105 


ios=0 

do  while  (ios.eq.O) 

read(30, 210, iostat=ios, end=1030)  aline 
210  format (a80) 

call  jplen (aline, i) 
write{40,*)  aline(l;i) 
end  do 

1030  continue 

close  (30) 
close  (40) 

C  Write  new  pafec  data  file  with  new  nodal  coordinates  added 
f ilename= 'JOB . DAT ' 
open  (unit=30, file=filename) 

open  (unit=40, file='TEMP.DAT' ,err=1040,status='old' ) 

ios=0 

do  while  (ios.eq.O) 

read  (40, 210, iostat=ios, end=1050)  aline 
call  jplen (aline, linelen) 
call  pulstrip (aline, linelen) 
if  (aline(l:8) .eq. 'PAFBLOCK')  then 
do  i=l, nodes 

write (30,225)  n (i,  j )  ,  x (i, j ) , y (i, j ) 

225  format (i3, 2fl0 . 3) 

end  do 

write(30,215)  'C 
215  format (al) 
end  if 

write ( 30 , 210 )  aline 
end  do 

1050  continue 
close (30) 
close (40) 

C  Increment  iteration  number  and  write  to  file 
j=j+l 

open (unit=10 , f ile= ' ITNUM.DAT ' ) 
write (10,235)  j 
close  ( 10 ) 

stop 

1010  write (6,*)  'Could  not  open  PUPPIES  output  file' 
stop 

1020  write (6,*)  'Could  not  open  Pafec  data  file' 
stop 

1040  write (6,*)  'Could  not  open  TEMP.DAT' 
end 

The  following  subroutines  by  courtesy  of  J.  Paul  AED,  AMRL. 
SUBROUTINE  JPLEN ( STRG , I LEN ) 

C  THIS  SUBROUTINE  OBTAINS  THE  LENGTH  OF  A  STRING 

C  INPUT:  STRG  -  STRING  TO  HAVE  LENGTH  DETERMINED 

C  OUTPUT:  ILEN  -  LENGTH  OF  STRING 

C 

CHARACTER* ( * )  STRG 
C 

C - PIND  LENGTH  OF  STRING. 

ILEN  =  LEN (STRG) 

10  IF  (STRG (ILEN: ILEN) .EQ. '  ')  THEN 

ILEN  =  ILEN  -  1 
IF  (ILEN.GT.O)  GOTO  10 
END  IF 
C 

c _ CHECK  FOR  FIRST  NULL  CHARACTER  IN  CASE  THERE  ARE  ONES 

DO  L1=1,ILEN 

N  =  ICHAR(STRG(L1:L1) ) 
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IF  (N.EQ.O)  THEN 
ILEN  =  LI  -  1 
GOTO  20 
END  IF 
END  DO 

20  CONTINUE 

C 

RETURN 

END 

SUBROUTINE  PULSTRIP (STRG, ILEN) 

C  THIS  SUBROUTINE  STRIPS  LEADING  BLANK  SPACES  FROM  A 

C  STRING 

C  INPUT  STRG  -  CHARACTER  STRING  BE  STRIPPED 

C  OUTPUT  ILEN  -  NUMBER  OF  CHARACTERS  IN  STRING 

CHARACTER* (* )  STRG,  CTEMP*300 
C 

CTEMP  =  STRG (1: ILEN) 

C 

Q - LOOP  TO  DETERMINE  FIRST  NON  BLANK  CHARACTER 

ILOC  =  1 

10  IF  (CTEMP (ILOC: ILOC) .EQ. '  ')  THEN 

ILOC  =  ILOC  +  1 
IF  ( ILOC. LT. ILEN)  GOTO  10 
END  IF 
C 

0 - RE-CONSTRUCT  STRING  IF  BLANK  CHARACTERS  EXIST 

IF  (ILOC.GT.l)  THEN 

STRG  =  CTEMP (ILOC: ILEN) 

ILEN  =  ILEN  -  ILOC  +  1 
END  IF 
C 

RETURN 

END 


9.4  File  names  and  descriptions 

JOB.DAT  This  is  the  name  of  the  PAFEC  data  file.  A  new  set  of  boimdary  nodes 
are  added  to  the  nodes  module  at  the  end  of  each  iteration.  Other  files  of  the  form 
JOB.*  are  created  by  PAFEC  in  the  process  of  solving  the  model  and  deleted  by  the 
PAFDEL  routine  at  the  start  of  the  next  iteration.  It  should  be  noted  that  PAFDEL  is  a 
pre-existing  in-house  routine  for  deleting  PAFEC  output  files  in  preparation  for  a  new 
PAFEC  run. 

JOBi.BS  These  are  the  PAFEC  binary  storage  files  for  the  model  including  loads 
and  restraints,  where  i  is  the  iteration  number.  For  each  iteration  they  are  copied  from 
the  active  JOB.BS.  They  do  not  contain  stress  results.  They  are  saved  for  later  plotting 
the  shape  at  any  iteration. 

JOBOPT.LOG  Stresses,  coordinates  and  other  key  quantities  for  each  iteration  are 
stored  in  this  file.  It  is  the  best  record  of  the  progress  of  a  solution. 

PUP.JNL  This  is  a  journal  file  for  a  data  extraction  program  called  PUPPIES.  It 
contains  the  list  of  node  numbers  that  are  to  be  moved  as  well  as  commands  to  extract 
average  nodal  principal  stresses. 
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JOB.PUPPIES  When  PUPPIES  is  executed  using  the  journal  file  PUP.JNL  a  list  of 
nodal  coordinates  and  stresses  is  written  to  this  file.  This  file  is  also  removed  prior  to 
the  next  iteration.  It  should  be  noted  that  PUPPIES  is  a  pre-existing  in-house  routine 
for  extracting  particular  PAFEC  output  in  numeric  form. 

ITNUM.DAT  Contains  the  iteration  number  in  i2  format. 

ITNUM.ONE  Is  copied  to  ITNUM.DAT  to  initialise  the  iteration  number.  Contains  1 
in  same  format  as  above. 

SPARE.DAT  The  initial  PAFEC  data  file  is  stored  here.  It  is  copied  to  JOB.DAT  to  re- 
mitiaHse  the  model  before  restarting. 

JOB.SOLVED  This  file  is  created  by  the  FORTRAN  program  when  convergence  is 
satisfied.  The  shell  script  terminates  when  this  file  is  present  in  the  directory.  The  file 
never  contains  any  data. 

TEMP.DAT  Used  for  temporary  storage  of  the  PAFEC  data  file  JOB.DAT. 


9.5  Shell  script  to  prepare  files  for  restarting  an  optimisation 

rm  JOB* . * 
rm  * . LOG 
rm  TEMP. * 

cp  SPARE.DAT  JOB.DAT 
cp  ITNUM.ONE  ITNUM.DAT 
F77  mn  itin.  f 

mn.f  is  the  fortran  program  that  moves  the  nodes  and  is  recompiled  between  runs  to 
incorporate  any  changes,  e.g  step  size. 
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